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1.0  INTRODUCTION 


In  a  previous  report^  it  was  shown  that  a  code  based  on  the  bihamonic 
formulation  of  the  Navier-Stokes  equations  with  the  corresponding  Newtonlzed 
difference  equations  solved  by  a  direct  band  solver  could  be  quite  competi¬ 
tive  in  situations  where  considerably  more  grid  points  are  used  in  one 
direction  than  in  the  other,  such  as  flow  in  a  long  channel.  The  main 
advantages  of  this  approach  are  that  it  is  straightforward  and  robust.  Both 
of  these  terms  are  subjective  of  course. 

To  justify  the  first  we  may  summarize  the  approach  as  follows:  From 
derivative  coefficients  generated  for  the  x  and  y  variable  grids  we  can 
immediately  write  down  the  difference  equations  |s*etTon  3)  and  it  is  then  a 
simple  matter  to  write  down  the  Newton  equations.  The  coefficients  are  stored 
by  diagonals  as  required  by  the  band  solver  and  the  righthand  sides  computed, 
then  the  Newton  corrections  are  obtained  by  calls  to  the  band  solver.  Also 
the  uniformity  of  the  approach  used  for  introducing  the  boundary  conditions  in 
section  4  simplifies  the  treatment  of  a  wide  variety  of  boundary  conditions. 

A- 

By  robustness  we  mean  less  likely  to  fail  in  a  wide  variety  of  situations 
where  other  methods,  possibly  more  efficient  but  requiring  careful  tuning, 
may  be  unable  to  produce  a  solution  at  all.  In  a  situation  where  one  has  an 
efficient  well  tuned  method  being  applied  to  a  familiar  problem,  it  is  dif¬ 
ficult  to  tolerate  the  large  storage  requirements  and  possibly  longer  computer 
times  needed  by  this  approach.  But  In  situations  where  an  unfamiliar  and 
awkward  problem  has  to  be  solved  quickly  without  too  much  tuning,  robustness 
can  be  a  very  welcome  property.  We  are  mainly  thi eking  of  the  avoidance  of 
the  all  Reynolds  number  problem,  which  was  demonstrated  in  reference  1.  This 
has  advantages  when  a  solution  of  some  sort  is  required  for  very  large  Reynolds 

numbers,  and  also  when  a  solution  is  required  for  a  moderate  Reynolds  number 

2 

on  a  coarse  grid,  for  example,  when  using  h  -extrapolation. 

As  discussed  in  reference  1,  the  main  drawback  of  the  approach  Is  the 
magnitude  of  the  storage  requirements.  Although  this  could  be  overcome  by 
making  use  of  direct-access  storage  devices.  It  was  felt  that  other  means 
of  reducing  the  storage  requirements  should  be  explored  first,  especially  if 
they  may  lead  to  a  reduction  In  computing  time  as  well.  The  present  report 
documents  progress  on  this  in  four  main  directions.  Also  in  two  appendices 
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tentative  work  is  reported  which  carried  out  an  alternative  to  the  bi harmonic 
scheme. 


The  bi harmonic  code  described  in  reference  1  was  restricted  to  laminar 
flows  and  uniform  grids.  The  code  for  which  improvements  are  being  considered 
is  a  more  general  version  of  that  laminar  code  in  which  turbulence  can  be 
taken  into  account  by  using  an  eddy  viscosity.  The  differential  equations 
for  this  form  of  the  stream  function  equation  are  derived  in  section  2. 

With  regard  to  the  improvements  to  this  code,  we  first  consider  the 
introduction  of  variable  grids  in  the  x  and  y  directions.  For  this  we  need 
formulas  for  higher  derivatives  up  to  the  fourth  in  which  the  grid  spacing  is 
not  assumed  to  be  uniform.  These  are  derived  in  section  3.  We  also  derive 
the  modifications  to  the  coefficients  in  the  Newton  matrix  that  this  general¬ 
ization  will  engender  in  section  5.  This  improvement  is  fully  implemented  and 
tested,  and  has  been  applied  to  the  test  problems  considered  in  sections  8  and 
9. 


The  second  technique  that  we  explore  is  to  tailor  the  band  solver  algo¬ 
rithm  so  as  to  take  full  advantage  of  the  considerable  number  of  zeros  that 
occur  within  the  band  of  the  matrix.  The  most  obvious  way  of  doing  this  is 
to  perform  the  LU  decomposition  by  a  block-elimination  scheme.  This  is 
considered  in  section  6,  where  the  formulas  for  the  forward  and  backward 
sweeps  are  derived.  Experiments  with  related  ADI  schemes  are  described  in 
Appendices  A  and  B. 

A  third  technique  that  could  be  used  to  introduce  some  flexibility  as  well 
as  to  reduce  storage  requirements  and/or  computer  time  is  to  divide  the  flow 
region  into  subregions  and  use  the  direct  LU  decomposition  procedure  on  each  of 
the  smaller  subregions.  Storage  requirements  for  the  decomposed  matrix  for 
the  smaller  subregions  will  be  much  less  at  the  expense  of  an  iterative  pro¬ 
cedure  to  couple  the  regions  together.  If  the  number  of  subregions  Is  fairly 
small,  this  outer  Iteration  should  converge  quickly  and,  since  the  Inner  iter¬ 
ations  for  the  subregions  should  be  much  faster.  It  may  be  possible  to  reduce 
computer  time  as  well  If  the  accuracy  aimed  at  Is  not  too  great.  A  simple 
scheme  for  linking  regions  Is  described  In  section  7,  but  there  is  a  potential 
here  for  more  sophisticated  schemes  In  which  one  could  allow  different  grid 
structures  in  different  regions,  thus  allowing  one  to  take  account  of  flow 
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structure  in  a  more  realistic  way.  This  could  be  particularly  advantageous 
for  large  Reynolds  numbers.  The  Implementation  of  a  simple  version  of  this 
technique  is  currently  being  considered  in  order  to  test  the  effectiveness 
of  the  basic  philosophy. 

These  first  three  techniques  are  quite  general  in  that  they  can  be  used 
for  whatever  problem  we  are  considering.  We  now  restrict  attention  to  channel 
flows  and,  in  particular,  the  sudden  expansion  configuration.  For  large 
Reynolds  numbers,  we  may  neglect  the  upstream  influence  in  the  entrance  channel 
and  assume  that  the  flow  is  fully  developed  right  up  to  the  entrance.  Thus  we 
may  take  the  entrance  conditions  In  the  expanded  channel  as  a  parabolic  profile 
over  the  middle  section  and  no  slip  over  the  remainder.  We  may  accommodate 
the  rapid  changes  from  these  upstream  conditions  by  taking  a  fine  x  grid  close 
to  the  entrance.  We  may  also  choose  a  fine  y  grid  near  the  walls  to  take 
account  of  the  boundary  layers  there.  This  means  that  we  will  have  a  fine  grid 
in  both  directions  in  the  corners,  which  should  take  care  of  any  problems  in 
the  corners.  The  chief  remaining  problem  area  is  the  slow  approach  to  fully 
developed  conditions  far  downstream.  This  can  be  dealt  with  to  some  extent 
by  choosing  an  expanding  x  grid  at  large  distances  from  the  entrance  region. 
However,  the  analysis  of  perturbations  from  fully  developed  conditions  due  to 
Wilson  shows  that  the  decay,  although  exponential,  can  be  very  slow  for  large 
Reynolds  number.  Thus  it  seems  very  worthwhile  to  make  use  of  the  asymptotic 
relations  derived  by  Wilson.  This  Is  considered  in  section  8.  We  also 
compare  its  effectiveness  with  a  simpler  downstream  condition,  which  we  call 
the  boundary-layer  condition.  This  may  be  as  effective  for  large  Reynolds 
numbers  and  would  be  especially  useful  for  turbulent  flows  since  we  do  not 
have  to  know  the  fully-developed  profile.  Results  for  turbulent,  two- 
dimensional  plane  channel  flow  are  presented  and  discussed  in  section  9. 


2.0  EQUATION  FOR  STREAM  FUNCTION.  INCORPORATION  OF  TURBULENCE  MODEL. 
The  turbulent  Navier-Stokes  equations  read 

ui!U  v|u  =_l|E  +  lIlxx  +l!!sn_i_(72)  _3_  (jjTyrv 
ax  ay  p  3x  p  ax  p  ay  ax  '  '  ay  '  ' 


uiv+  av  =  _la£  +  l^  +  l!^_a_  (-rV-)  _  a_ 
u  ax  ay  p  ay  p  ax  p  ay  ax  '  '  ay  1  ' 


(2.1) 


Since  txx  =  2p(au/ax),  t  =  u(au/ay)  yfav/ax),  Tyy  =  2P(av/ay),  these 
can  be  written 


u  Ju  +  v  =  _  1  3R  +  JL  lii _ -^7]  +  a_  ,  {au  +  av} 
ax  ay  p  ax  ax  1  ax  u  J  ay  lv'ay  ax'  u  1 


uiv  +  viv  =  _11E  +  1_  [v(lJUiV)  -XTVr-]  +  [2v  » 1-7^] 

ax  ay  p  ay  ax  1  vay  ax'  '  ay  1  ay  J 


(2.2) 


Making  the  assumption  that  the  Reynolds  stresses  can  be  modelled  by  an  eddy- 
viscosity  approach,  we  write 

2v  a . ^ *  g  ?x.  ,<£ ♦  *>  - .  (v  ♦  «„><£  *  *) 


2vM_^.2(vt«m)  B, 


Then  if 


nr' ay  ax' 

(2.3) 

(2.4) 


b  -  0  +  «J)»  where  4  -  em/v  (2.4) 

and  we  Introduce  nondimensional  variables,  we  obtain 

uiyi  +  v  —  =  — -^  +  i—  (2b  — )  +  i  —  [b(—  +  ■^■)1 

u  ax  ay  ax  E  ax  '  ax'  1?  ay  1  'ay  ax'J 

(2.5) 

u  3v  av  B  _  a£+  1  a_  [b(au  +  av}]  +  1  a_  (2b  |v} 

u  ax  ay  ay  R  ax  1  vay  ax"  R  ay  '  ay 

Eliminating  the  pressure,  we  obtain 

„  a  ,au  av%  .  „  a  rau  av ,  „  1  a2  Lh/j>u  _  avJ  +  /a2  _  a2  Jhfiil  +  l£)]l 
u  a?  a?}  v  ay  ay  ~  ax  Tf  ( axay  [Zb{n  W\  ^  ^  3X  i| 
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which  in  terms  of  the  nondimensional  stream  function  f  reads 


<>1*2  3f 
ay  ax 


af  „2  af  _  1 
ax  v  ay  ~  £ 


a2 

4b  *2f 

axay 

axay 

(a2 

ay 


a2  ) 

I?’ 


bfa2f 

b(~T 

ay 


a2f\ 

I? 


Differentiating  out  we  obtain 


(2.6) 


4(bf  )  =  4bf  +  4b  f  +  4b  f  +  4b  f 

xy  xy  xxyy  x  xyy  y  xxy  xy  xy 


[b^fyy  fxx^yy  b  (fyy  *Wyy  +  2Vfyy  fxx^y  4  byy^fyy  fxx^ 
'[b^fyy  ~  fxx^xx  “  “b^fyy  _  fxx^xx  -  2bx^‘fyy  +  fxx^x  “  bxx^fyy  ~  fxx^ 


fyV  fx  V  fy  £  [b^fxxxx  +  2fxxyy  +  fyyyy^  +  2bxv  fx  +  2byv  fy  +  4bxyfxy 

+  (byy-bxxHfyy-fxx^ 


Thus  finally  the  "turbulent"  stream- function  equation  reads 

V2fx  -  f/fy  ■  l  *  "A*  +  (byy  -  bxx>‘fyy  -  fxx> 

*  4bxyV  <2-7> 
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3.0  DIFFERENCE  FORMULATION  WITH  VARIABLE  GRID 


We  set  up  a  rectangular  grid  covering  the  basic  rectangle  x^  <  x  <  ^ 

y  •„  <  y  <  y„,u  with  M  grid  lines  in  the  x-di  recti  on  and  N  in  the 
mi  n  max 

y-direction.  Actually  we  also  store  two  exterior  x  and  y  grid  values  at 
each  end  because  we  need  them  in  computing  coefficients  for  derivatives 
actually  on  the  boundaries.  The  grid  spacing  in  each  direction  may  be 
non uni  form. 


Since  the  stream  function  equation  involves  drivatives  of  up  to  the  fourth 
order,  we  require  formulas  for  higher  derivatives  at  the  grid  points  of 
a  grid  with  nonuniform  spacing.  We  restrict  attention  to  the  x  grid 
and  indicate  only  briefly  the  modifications  in  notation  needed  to  transform 
to  the  y  grid.  We  can  represent  to  second-order  accuracy  the  1st  and  2nd 
derivatives  in  terms  of  three  neighboring  function  values  by  locally  fitting 
a  quadratic.  Thus  we  find  that 


rvvi  i 

'hi 

-hi-ii 

hi-l/h1 

hi+hi-i  _ 

fi-i + 

*hi-i. 

f1 

h1+h1-l _ 

(3.1) 


p/h1-l  1 

f 

2 

f  + 

'  2/h1  ‘ 

*  hi'hi-l 

[v>v7] 

ri-l 

T1  + 

h1+hi-l. 

f1+l  '  * — 3 — 

where  h^  *  x^-x^.  From  the  error  terms  we  see  that  the  first  formula  Is 
clearly  second  order  and  so  Is  the  second  If  the  rate  of  change  of  grid  spacing 
with  x  is  of  order  unity. 

For  the  3rd  and  4th  derivatives  we  need  to  fit  a  quartlc  locally,  so  these 
formulas  will  Involve  five  neighboring  function  values  ^i_i»  ^i+i» 

f^+2«  It  will  be  convenient  to  adopt  Immediately  a  notation  that  can  be 
transcribed  directly  into  the  Fortran  code,  so  we  write 

f.|’  *  a2^  )^1-1  +  a3^  )^i  +  a4^ )^i+i 

*  c2(i)f1-l  +  c3(i)f1  +  c4(i)f1+l  (3.2) 

ff  -  P^Off.2  +  P2(1)fi_i  +  P3(Dfi  +  P4(1)fi+1  ♦  P5(1)f1+2 

T  '  rl ^ >fl-2  +  r2^)f1-l  +  '3(1>f  1  +  r4^1>fi+l  +  r5^)f1+2 
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where  the  coefficient  notation  b  k ( i ) ,  ( i ) »  q  k ( i ) ,  sk(i)  has  been  reserved 

for  the  corresponding  coefficients  for  derivatives  in  the  y  direction. 

Formulas  for  the  a^l)  and  c ^ ( i )  have  been  given  above  and  can  easily 
be  derived  from  Taylor  expansions.  However,  for  the  higher  derivatives  the 
algebra  is  more  complicated  and  it  Is  convenient  to  work  more  systematically  in 
terms  of  divided  differences.  Thus  Newton's  Interpolating  quartic  reads 


f(x)  =  A0  +  A'j  ( X-Xj  )  +  a2(x-x1)(x-x2)  +  ...  +  aa[x-x,)...(x-xa)  (3.3) 


4V*-*'! , 


where  we  have  for  convenience  taken  i  =  3;  a  shift  to  the  general  local  grid 

point  can  easily  be  made  whenever  required.  The  coefficients  occurring  here 

are  the  divided  differences  [2] 

net  n+1  -| 

%■  H  J  <*j-xk>"1  ,(xj>  (3-4) 

J=1  ^  k« 

To  derive  the  higher  derivative  formulas  we  write 

9  3  2 

f(x)  *  AQ  +  A-j(x-Xj)  +  A2[x4  -  (Xt+Xj^X  +  X.,X9]  +  A,[x  -  (x^Xp+X^Jx  +  . 


rA2'A  A1A  2J  3L 


^4[x  ~  (xi+x2+x3+x4)x  +  . ..] 


f'"(x)  =  6a3  +  6a4(4x-x1-x2-x3-x4),  f""(x)  =  24a4, 

and  for  the  values  at  the  "central"  grid  point  x3  we  obtain 
f'"(x3)  =  6[a3  +  a4(3x3-x1-x2-x4)] 
f""(x3)  -  24a4 

Thus  the  coefficients  r. (3)  are  Immediate  from  the  formulas  for  A  : 

5  K  n 

rt(3)  =24  TT  (x^-x,)'1 


(3.5) 


(3.6) 


1 
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j'v 


and  the  p. ( 3 )  can  be  obtained  from  them  by 

K  4 

Pk(3)  =  6  TF  (Vx  )-1  +  ^3x3-x-|-x2-x4)rk(3)  for  k  =  1,  2,  3,  4 
J  1 
ji*k 

and  finally 

Pg(3)  =  (3x3-x^-X2"X^)rg(3) . 

These  formulas  for  the  Pk(3)  can  be  summarized  more  conveniently  as 

pk(3)  =  l  (7+  xk)rk(3),  k  =  1,  2 . 5,  (3.7) 

where 

x  =  3x3  —  x-j  —  X2  —  x4  —  Xg.  (3.8) 

These  coefficients  and  the  corresponding  ones  for  the  y  grid  are  assumed 
to  be  pre-calculated  and  stored  for  each  grid  line.  To  obtain  the  difference 
equations  for  the  grid  stream  function  F.  .  we  first  write  down  the  formulas 
centered  on  the  grid  point  (i,j)  for  the  various  derivatives  and  cross 
derivatives  that  occur,  omitting  for  simplicity  the  i  or  j  index  associ¬ 
ated  with  the  derivative  coefficients  —  for  ak,  ck>  pk>  r^  this  will  always 


be  i 

and  for 

bk 

,  dk. 

V 

,  it  will  always  be 

j: 

f 

xxxx 

=  rlFi-2 

.  + 
.J 

r2Fi 

-1  J 

+ 

r3Fi,j  +  r4Fi+l  ,0  * 

r5Fi+2,j 

Fxxyy 

=  c2d2Fi 

■1  J- 

-1  + 

c2d3l 

F. 

1- 

-1  J  +  C2d4Fi  - 1  J+l  +  +  c4d4Fi+l 

J+l 

fyyyy 

=  SlFi,j- 

-2  + 

s2Fi 

J-l 

+ 

S3Fi  ,j  +  s4F1  ,j+l  + 

S5Fi  J+2 

f 

XXX 

=  PlFi-2 

.  + 
,J 

P2Fi 

-1  J 

+ 

P3F1,<  *  p4Fi+7  ,j  * 

P5Fi+2  J 

f 

xxy 

=  c2b2Fi_ 

■1  J- 

•1  + 

c2b3l 

F. 

1- 

■1J  +  C2b4Fi-l  J+l  +  +  c4b4Fi+l 

J+l 

fxyy 

=  a£d2Fi- 

■l»j- 

■1  + 

a2d3* 

i- 

■l,j  +  a2d4Fi-l  ,j+l  4 

'  •”  +  a4d4Fi+l 

J+l 

fyyy 

=  qlFi,j- 

■2  + 

q2Fi 

J-l 

+ 

fl3Fi,j  4  "4Ff.j+l  + 

q5Fi ,j+2 
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fxx  =  c2Fi-l,j  +  c3Fi,j  +  C4Fi+l ,j 
fyy  =  d2Fi,j-1  +  d3Fi,j  +  d4Fi  ,j+l 


Fx  a2Fi-l,j  +  a3Fi,j  +  a4Fi+l,j 


f  =  b0F.  .  •>  +  b.F.  .  +  b.F.  .  , 
y  2  i  ,j-l  3  i  ,J  4  i  ,j+l 


Corresponding  formulas  can  be  written  down  for  the  derivatives  of  the  eddy 

viscosity  b  which  occur  in  the  stream  function  equation,  namely  b  ,  b  , 

xx  yy 

b  ,  b  . 
x  y 


If  we  divide  through  by  the  eddy  viscosity,  the  stream  function  equation 


reads 


4f  +  2bx  Rfy  y2f  +  2by  +  Rfx  2.  +  4  f  +  bxx  b,y 
E>  x  b  y  b  xy  E~* 


<fxx  -  fyy)  ■ 


which  can  be  written  for  convenience  in  deriving  the  Newton  equations  as 


E  =  0 


(3.10) 


where 


E  fxxxx  +  2fxxyy  +  fyyyy  +  A^fxxx  +  fxyy;  +  B^fxxy  +  fyyy^  +  Cfxy 

+  D<fxx  -  V 


(3.11) 


2b..  -  Rf. 


,  B  = 


2b..  +  Rf 


,  C  =  — ,  D  = 


b  -  b 
xx  y 


(3.12) 


The  nonlinear  difference  equations  are  obtained  by  replacing  the  f  derivatives 
here  by  the  variable  grid  formulas  derived  above  for  each  interior  grid  point 
and  adding  boundary  condition  equations;  see  next  section. 


4.0  BOUNDARY  CONDITIONS 


We  need  to  specify  two  boundary  conditions  on  each  of  the  four  sides. 

These  usually  involve  f  and  its  normal  derivatives  f  or  f  and  may 

n  nn 

be  linear  or  nonlinear.  Examples  of  typical  linear  boundary  conditions  are  the 

following.  Along  an  inlet  boundary  we  usually  specify  u  and  v,  i.e.  f 

and  its  normal  derivative.  We  also  specify  u  and  v  on  a  fixed  wall  when 

the  flow  is  laminat .  Along  a  line  of  symmetry  we  can  specify  f  constant  and 

zero  second  normal  derivative.  On  a  downstream  boundary  several  pairs  of 

boundary  conditions  can  be  used,  which  represent  the  approach  to  fully 

developed  conditions  with  varying  degrees  of  effectiveness.  Examples  of  linear 

ones  are  fn  =  f  =  0,  which  essentially  says  that  f  does  not  vary  with 

the  normal  coordinate  n,  and  f  =  f  =0,  which  we  call  the  boundary- 

nn  nnn 

layer  condition.  A  more  sophisticated  downstream  condition  is  that  the  approach 
to  the  downstream  profile  is  exponential  (see  Wilson,  1969,  for  the  laminar 
case).  This  leads  to  an  example  of  a  nonlinear  condition  if  we  do  not  know 
the  decay  rate,  as  discussed  below.  Another  example  of  a  nonlinear  condition  is 
the  law-of-the-wall  condition,  which  may  be  used  to  replace  the  fp  =  0  con¬ 
dition  on  a  fixed  wall.  This  will  be  discussed  further  in  section  9  where 
specific  turbulent  problems  are  treated. 


The  main  point  we  wish  to  make  about  the  above  list  of  pairs  of  boundary 
conditions  is  that  they  can  all  be  represented  in  the  interests  of  uniformity 
as  a  pair  of  the  form 

f0  =  q(frf2).  fT  «  r(f2,f3)  (4.1) 

where  the  subscripts  count  grid  lines  inwards  across  the  boundary  with  the 
boundary  itself  characterized  by  the  subscript  1.  Thus  f  denotes  the 
fictitious  exterior  value  which  is  not  actually  stored.  The  functions  q  and 
r  may  be  linear  or  nonlinear.  In  either  case  their  partial  derivatives,  which 
will  be  needed  in  the  Newton  linearization,  will  be  denoted  by 


(4.2) 


Thus,  for  example,  the  linearized  version  of  the  second  condition  would  read 


6 f i  -  c«f2  -  d«f3  =  r(f2,f3)  -  f1 


(4.3) 


10 


This  form  is,  in  fact,  used  directly  in  the  next  section,  whereas  the  first 
condition  is  used  essentially  to  substitute  for  the  unstored  exterior  values. 


For  incorporation  in  the  system  of  algebraic  equations,  these  conditions 

have  to  be  particularized  for  each  of  the  four  sides,  which  we  do  by  attaching 

a  superscript  N,  S,  E,  W  corresponding  to  the  North,  South,  East  or  West 

side.  Also  they  have  to  be  written  in  terms  of  the  local  grid  stream  function 

values  F.  -,  where  i  =  1  or  M  along  the  South  or  North  side  and  j  =  1 

or  N  along  the  West  or  East  side,  respectively.  Since  the  q,  a,  b,  r,  c,  d 

may  vary  along  the  sides,  they  are  all,  in  the  interests  of  uniformity,  stored 

as  one-dimensional  arrays.  As  an  example,  suppose  the  downstream  condition 

f  =  f _  =  0  is  to  be  imposed  on  the  North  side.  The  first  condition  f„„  =  0 

nn  nnn  r  nn 

can  be  represented  directly  by  using  the  variable  grid  derivative  coefficients 
corresponding  to  the  grid  line  j  =  N.  The  second  condition  says  essentially 
that  the  second  derivative  is  not  changing,  so  we  can  represent  it  by  imposing 
f  =0  also  on  the  first  interior  grid  line  j  =  N-l.  The  pair  of  conditions 
therefore  read 


¥N>F<,N-1  *  d3<N>F1.N  +  d4<N»Fi.H+l  -  0 
d2(N-l)F1iN.2,d3(B-l)F1iK.1*d4(H-l)Fj>N=0 
or  in  the  notation  of  (4.1) 

Fi  ,N+1  =  qN(Fi,N,Fi,N-l^  Fi  ,N  =  r^Fi  ,N-1  ,Fi  ,N-2^ 


where 


with 


^F1,N-W  =  aVt<  +  bN,:l,N-l 

rNlFi,N-l'Ft,N-2'  °  cNp1,N-l  +  dNp1 ,N-2 


aN  =  -d3(N)/d4(N),  cN  =  -d3(N-l)/d4(N-l) 

bN  =  -d2(N)/d4(N),  dN  =  -d2(N-l)/d4(N-l) 

Details  for  the  other  conditions  will  be  given  where  they  actually  arise. 


(4.4) 


(4.5) 


(4.6) 


(4.7) 
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5.0  NEWTON  MATRIX  AND  SOLUTION  BY  BAND  SOLVER 


Let  F  be  the  vector  of  MN  unknown  grid  stream  function  values  F.., 

*  J 

i  =  1,  ....  M,  j  =  1,  ....  N,  i.e.  the  interior  values  plus  the  boundary 
values,  but  not  Including  the  exterior  values.  We  obtain  MN  equations  for 
these  unknowns  by  writing  down  the  grid  stream  function  equation  for  all  the 
interior  grid  points  (MN  -  2M  —  2N  +  4  equations)  with  exterior  values 
assumed  to  be  expressed  in  terms  of  the  F. .  by  means  of  the  first  boundary 
condition  along  each  side,  adding  2(M  -  2)  +  2(N  -  2)  equations  representing 
the  second  boundary  condition  on  each  side  and  finally  adding  four  equations 
for  the  corners.  The  last  are  obtained  by  selecting  the  second  boundary  con¬ 
dition  from  either  of  the  two  sides  meeting  at  the  corner.  In  our  code  we 
optionally  choose  these  sides  to  be  the  North  and  South  sides. 

Let  6j.(F),  i  =  1,  ...»  M,  j  »  1,  ...»  N  be  functions  defined  so  that 

1 J 

the  equations 

3id  (F)  =  0,  1-1 . M,  j  =  1 . N  (5.1) 

represent  the  above  MN  equations.  Then  with  Newtons  method,  if  F  does 
not  satisfy  (5.1)  accurately  enough,  it  is  corrected  to  F  +  $  where  ^ 
satisfies 


“u  ‘N+y.j+v 

I  v»  | + 1  v  |  <2 

the  being  the  partial  derivatives  of  with  respect  to  the  that 

it  depends  on.  From  (3.9),  (3.11),  (3.12)  we  find  that  with  the  abbreviations 


9x  R^xxx  +  fxyy^» 


9,  '  R(fxxy  +  W/b 


we  can  write  down  immediately 


a_2*  f*i  +  Ap.j 
a2  :  r5  +  ^5 

°o2:  S1  +  Bql 
2 

a*;  :  sc  +  Bac 
0  5  5 
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a_|:  a2(Ad2  +  Cb2)  +  c2(2d2  +  Bd2 ) 

a^:  a4(Ad2  +  Cb2)  +  c4(2d2  +  Bd2) 

o)l :  a2(Ad4  +  ^4)  +  c2(2d4  +  ^4) 

«]'•  a4(Ad4  +  ^d4 )  +  c4(2d4  +  Bd4) 

«?,:  a2(Ad3  +  Cb3  +  g  )  +  c2(2d3  +  Bd3  +  D)  +  r?  +  Ap2  (5.2) 

a°  :  a4(Ad3  +  Cb3  +  gy)  +  c4(2d3  +  Bd3  +  D)  +  r4  +  Ap4 

o0^:  b2(Bc3  +  Ca3  "  9X)  +  d2^2c3  +  Aa3  -  D)  +  s2  +  Bq2 

°o  :  b4^Bc3  +  Ca3  ~  9X}  +  d4^2c3  +  Aa:  “  +  s4  +  Bq4 

“o  :  a3^Ad3  +  Cb3  +  9y)  +  c3^2d3  +  Bb3  +  0)  +  r3  +  Ap3 

-  d3D  -  b3gx  +  S3  +  Bq3 


When  i  =  2  or  M  -  1  or  j  =  2  or  N  —  1 ,  the  a2,  a2,  a”2  or 
a°2  would  multiply  4> * s  corresponding  to  exterior  values,  so  modifications 
must  be  made  to  take  into  account  the  relevant  first  boundary  conditions. 

For  example,  if  i  =  M-l  and  3  <_  j  <^N-2,  as  in  Fig.  1  the  relevant  first 
boundary  condition  would  be  of  the  form 

FM+1 , j  =  qE(FM,j’  FM- 1 , j  ^ 


with 


n,j 


The  derivatives  of  with  respect  to  FM  j  and  FM_u  then  becomes 
aj(i,j)  +  aE(j)a°(1,j)  and  a°(i,j)  +  bE(j)a°(1,j) 
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Figure  1.  A  typical  computational  molecule  at  a  point  where  an  exterior 
unstored  grid  value  Is  used.  The  dotted  lines  Indicate  the 
grid  values  connected  by  the  two  boundary  equations  (4.5). 

Stored  and  unstored  grid  values  in  the  calculation  are  Indicated 
by  •  and  o. 
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Thus  (*2  is  not  used  in  the  matrix,  but  is  used  to  modify  and  a°. 

Similar  modifications  must  be  made  near  the  other  boundaries.  For  points 
neighboring  both  boundaries,  two  exterior  values  are  involved  so  modif1ca- 
tions  must  be  made  for  each;  thus  a°  will  then  be  modified  twice. 

The  Newton  equations  for  the  second  boundary  conditions  are  added  as 
appropriate  and  have  the  form,  for  example,  of 

♦M,j  ~  cE*M-l,j  ~  dE*M-2,j  =  rE(FM-l,j‘  FM-2  ,j  )  "  FM,j 
when  i  =  M. 

There  are  various  ways  available  for  the  solution  of  the  linear  system 
formed  by  these  MN  equations.  Several  block  elimination  schemes  can  be 
formulated  and  also  iterative  schemes,  including  those  of  ADI  type.  If  one 
wishes  to  take  advantage  of  the  chord  method,  i.e.  the  simplified  Newton 
method  in  which  only  the  right-hand  sides  are  updated  at  each  iteration,  a 
scheme  in  which  an  LU  decomposition  is  performed  and  stored  would  seem  to  be 
preferable.  Block  elimination  schemes  can  be  designed  to  do  this,  but 
because  of  its  potentially  greater  stability,  espeically  for  large  Reynolds 
numbers,  we  chose  to  explore  first  the  practical  application  of  a  standard 
band  solver  with  partial  pivoting  for  stability  (this  may  be  especially 
important  if  we  work  entirely  in  IBM  single  precision). 

The  standard  way  of  organizing  the  equations  as  a  banded  system  is  to 

order  the  unknowns  by  rows,  i.e.  in  the  order  i »  •  •  • »  *frm  -| »  "fr]  2’  '  ’ ’ 

The  matrix  elements  for  the  system  (5.2)  then  has  the  appearance  shown  in 
Fig.  2.  To  avoid  confusion  the  suffices  attached  to  the  a’s  are  shown  on 

one  row  only,  and  the  modifications  to  the  a's  due  to  the  first  boundary 

conditions  are  only  indicated  by  position.  Thus  those  needing  modification  at 
the  West  and  East  boundaries  are  shown  with  a  bar,  those  at  the  North  and  South 
by  a  hat,  and  those  at  both  by  a  bar  and  a  hat. 

For  use  In  the  band  solver  the  diagonals  have  to  be  stored  as  the  columns 
of  a  rectangular  array  D,  say,  so  they  are  numbered  from  1  to  ^  =  4M+1. 

The  main  diagonals  of  the  other  Individual  blocks  also  have  salient  locations, 
which  we  therefore  denote  by  ^ ,  *2,  1^,  so  that  other  diagonals  can  be 
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referred  to  them  as  indicated  in  a  typical  row  in  Fig.  2.  An  equation  count  k 
is  also  shown  in  relation  to  the  successive  blocks,  together  with  the  right- 
hand  sides. 

Since  the  coefficients  are  typically  of  order  (grid  size)"^,  they 
are  much  larger  than  the  coefficients  in  the  boundary  equations.  The  solver 
therefore  interchanges  them  in  pivoting  for  size,  which  causes  loss  of 
accuracy.  Strictly  we  should  scale  all  the  equations  so  that  they  are  all 
fairly  uniform.  However,  we  find  that  the  natural  scaling  of  the  stream 
function  equations  is  adequate  provided  we  scale  up  the  boundary  equations 
sufficiently  so  that  they  dominate  and  are  not  interchanged. 

A  discussion  of  storage  estimates  and  operation  counts,  as  well  as 
iteration  strategies  was  given  in  reference  1  and  need  not  be  repeated  here, 
apart  from  specific  details  given  in  the  examples. 
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6.0  QUIN-BLOCK-DIAGONAL  SOLVER 


A  suite  of  subroutines  has  been  written  and  tested  for  the  LU  decompo¬ 
sition  and  solution  of  the  quin-block-diagonal  system  arising  in  the  Newton 
iteration  scheme.  The  system  has  the  form 
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where  the  blocks  are  m  x  m  square  matrices.  Here  we  have  chosen  the  safer 
form  of  decomposition  in  which  the  factors  are  exactly  triangular.  The 
diagonal  blocks  Cj  of  the  left-hand  factor  are  unit  lower  triangular. 

The  recurrence  relations  for  the  LU  decomposition  used  in  subroutine 
LUBL05  read 

1.  *•  Cj_2  =  A.  (3  <  j  <  n) 

2.  B'.  -  B.  -Aj  0\_2  (2  <  j  <  n) 

3-  cj  c"  ■  pj  <cj  - Bj  Dj-i  - #j  Ej.2'  "iii-i 

«•  S'  Dj  ■  pj(Dj  -  Bj  E5-i> 

5'  Cj  Ej  =  Pj  Ej  0  i  j  <  "-2) 
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wh-ire  the  Pj  are  permutation  matrices  introduced  by  the  interchanges  in 
the  LU  decomposition  of  C. . 

J 

A  subroutine  SP  has  been  written  for  subtracting  the  products  involved 
in  (2),  (3),  (4),  a  subroutine  TS  for  the  transposed  forward  substitutions 
in  (1)  and  (2),  and  a  subroutine  FS  for  the  forward  substitutions  in  (4)  and 
(5). 


A  subroutine  SOBLE  uses  the  decomposition  to  solve  the  block  equations. 
The  recurrence  relations  involved  in  the  foward  and  backward  substitutions 
read 

Forward:  =  P^,  =  P2(r2  -  ) 


Ciei 


Backward: 


Cu<b  =  e  , 
n^n  mi 


Cn-l£n-l 


7. 


clt 


J*j  °j  Dj®j+1  Ej®j+2 


(3  <  j  <  n) 

e  ,  -  D1  .e 
mi-1  n-Kn 

(1  <  J  <  n-2) 


A  vector  back  substitution  subroutine  BS  was  written  for  solving  the  equa¬ 
tion  in  the  backward  sweep.  Also  it  turned  out  to  be  convenient  to  write  a 
vector  subroutine  S2TV  for  subtracting  off  the  two  transformed  vectors 
which  occur  in  (6)  and  (7). 


The  above  subroutines  have  all  been  tested.  Some  significant  loss  in 
accuracy  was  observed  in  one  test,  but  this  is  believed  to  be  due  to  the 
particular  test  matrix  chosen  being  somewhat  ill-conditioned.  Moderate  loss 
in  accuracy  was  also  observed  in  another  test  where  the  matrix  was  of  blhar- 
monic  form,  so  It  seems  that  further  Investigation  Is  needed  before  the  use 
of  a  block  solver  Is  Implemented  in  the  present  code.  Also,  in  the  present 
form  the  saving  in  storage  Is  not  substantial.  Further  savings  could  be  made 
if  a  more  sophisticated  storage  scheme  was  adopted,  but  It  might  be  simpler 
to  try  the  alternative  pseudo-LU  decomposition. 
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7.0  ITERATION  BETWEEN  SUB-REGIONS 


Another  scheme  for  reducing  the  storage  requirements  Is  to  divide  the 
flow  region  Into  sub-regions  and  use  the  direct  LU  decomposition  procedure 
on  each  sub-region  in  turn.  Once  the  local  decomposition  has  been  performed. 

It  would  seem  worthwhile  to  iterate  the  simplified  Newton  (chord)  procedure 
several  times  before  proceeding  to  the  next  sub-region.  The  problem  of  connecting 
one  region  to  the  next  can  be  solved  very  simply  In  the  case  when  the  grid 
structure  is  the  same  for  all  sub-regions:  we  merely  use  the  band  matrix  for  the 
whole  region  to  provide  the  coefficients  for  the  connecting  relations  between 
values  on  either  side  of  the  dividing  line.  Of  course  the  remaining  band 
sub-matrices  provide  the  band  matrices  for  each  sub-region.  This  Is  Illustrated 
in  Figure  3  for  the  case  of  three  sub- regions.  The  scheme  would  apply  for 
any  very  long  band  matrix,  but  in  the  present  biharmonic  context  each  division 
would  normally  be  chosen  to  correspond  to  the  end  of  a  grid  line  and  the  overlap 
vectors  yk,  zk  would  then  correspond  to  the  two  neighboring  grid  lines  on  the 
other  side  of  the  boundary.  The  Iterative  procedure  we  would  use  would  be 
essentially  a  block  Gauss-Seldel ,  which  for  the  case  represented  In  the  diagram 
would  be  written 


AjX 

A2X 

A,x 


(1) 

1 

.  r,  -  0,4'-') 

(1) 

2 

-  r2-C2x<'-'>- 

(1) 

3 

=  r3  -  B3x^ 

,d> 


where  1  Is  the  Iteration  count  and  x^"^.  are  known  from  the 

previous  Iteration.  Note  that  because  of  the  structure  of  the  matrices  Bk>  Ck 
only  the  overlap  parts  of  the  vectors  xk*"^  are  Involved.  If  necessary, 
an  acceleration  parameter  could  be  Introduced  to  accelerate  convergence. 


Various  refinements  and  extensions  to  this  basic  scheme  can  be  envisaged. 
For  example.  If  the  sub- regions  are  chosen  so  that  within  each  the  number  of 
unknowns  In  one  direction  Is  more  than  In  the  other  they  should  be  ordered 
so  that  the  M  In  the  estimate  8M3N  for  the  LU  operation  count  should  be 
less  than  N.  A  generalization  that  would  require  rather  more  effort  to 
Implement  but  which  should  ultimately  be  very  effective  Is  to  reformulate  the 
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Figure  3.  Newton  matrix  for  case  of  three  subregions. 
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matching  conditions  at  the  sub-region  boundaries  in  such  a  way  that  different  grid 
structures  could  be  used  in  each  sub-region.  This  would  involve  interpolation 
and  perhaps  reformulation  in  terms  of  boundary  derivatives,  although  one  could 
probably  get  away  with  fairly  minor  modifications  to  the  present  variable  grid 
code  if  one  restricted  attention  to  the  situation  where  the  boundary  between 
regions  is  entirely  in  the  fluid,  i.e.  no  part  of  it  also  has  boundary  condi¬ 
tions  specified.  As  a  fairly  typical  example  we  may  consider  the  expanding 
channel  of  Figure  4,  where  the  flow  region  is  divided  into  three  sub-regions 
I,  II,  III.  The  simple  situation  envisaged  exists  for  the  boundary  between 
II  and  III  and  also  for  that  between  I  and  II  looked  at  from  Inside  I,  but 
not  from  inside  II. 

We  therefore  consider  first  solving  for  region  III  assuming  values  in  II 
are  known  from  a  previous  iteration.  We  imagine  the  Ill-grid  extended  into 
region  II  for  two  y-grid  lines.  The  x-grid  lines  do  not  match  so  we  interpolate 
to  obtain  known  values  on  the  two  grid  lines  outside  of  III.  Since  these  are 
known,  the  Newton  corrections  for  these  values  are  all  zero,  so  we  can  merely 
truncate  the  Newton  matrix  at  the  appropriate  point  and  not  complete  it  with 
boundary-condition  equations  and  modifications.  The  code  changes  required  for 
this  is  relatively  straiahtforward.  so  we  are  currently  considerina  implementing 
them  in  order  to  assess  the  effectiveness  of  the  basic  philosophy. 


Figure  4.  Expanding  channel  divided  into  subregions. 


It  is,  of  course,  entirely  feasible  to  deal  with  the  more  complicated 
situation,  but  since  more  sophisticated  code  changes  would  be  required,  it 
seems  reasonable  to  hold  this  extension  in  abeyance  for  the  immediate  future, 
pending  the  above  assessment. 
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8.0  LAMINAR  FLOW  APPLICATIONS 


We  first  tested  the  variable  grid  code  on  the  idealized  channel  flow  prob¬ 
lem  with  uniform  parallel  inlet  velocity  already  calculated  by  the  uniform  grid 
code1.  We  omit  details  of  the  application  of  the  code  since  only  minor  modifi¬ 
cations  were  needed  from  those  given  in  reference  1.  It  was  clear  that  the 
singularities  at  the  corner  could  be  dealt  with  more  satisfactorily  by  concen¬ 
trating  grid  points  near  them,  but  since  this  situation  is  somewhat  idealized, 
it  was  decided  not  to  put  too  much  effort  into  producing  accurate  results  over 
a  large  Reynolds  number  range  for  this  problem,  but  to  transfer  attention  to  the 
more  realistic  problem  of  the  sudden  expansion  channel.  Suffice  it  to  say  that 
some  experimentation  was  performed  on  choice  of  grid  structure  and  that  the 
original  choice  of  geometric  grid  with  step  h^  =  h^k1'1  (k  >  1)  was  ultimately 
abandoned  in  favor  of  a  quartic  grid  in  which  the  step  varied  cubical ly,  the 
maximum  and  minimum  of  the  cubic  being  at  the  boundaries,  the  center  of  the 
channel  and  the  wall  for  the  transverse  grid  and  at  the  last  station  and  the 
inlet  for  the  downstream  grid.  This  meant  that  for  the  transverse  grid  the  step 
became  large  but  fairly  uniform  far  downstream,  which  is  appropriate  for  the 
exponential  decay. 

We  turn  now  to  the  sudden  expansion  case  and  again  take  x  across  the 
channel  and  y  along  the  channel  as  in  Fig.  5  with  x  =  0  on  the  wall  and 
x  =  1/2  on  the  centerline.  We  assume  a  fully  developed  parabolic  profile  at 
the  inlet  over  a  central  opening  at  y  =  0.  This  starts  at  x  =  a,  where  we 
have  taken  a  =  1/3  but  it  can  be  varied.  The  choice  of  a  =  1/3  was  made 

4 

because  accurate  experimental  results  are  available  .  The  Reynolds  number 
quoted  for  these  results  was  Re  =  56  and  was  based  on  the  maximum  velocity 
of  a  profile  slightly  upstream  of  the  inlet.  Assuming  that  this  was  a  fully- 
developed  profile,  the  corresponding  Reynolds  number  in  our  notation  would  be 
R  =  41.3.  The  boundary  conditions  on  y  -  0  are 

f  ■  F„(x).  fy  -  0 

where 

0  If  0  <  x  <  a 

Fo(x)  =  (3  -jprfMTHH)2  1*  a<x<l/2 


24 


f  =  f  =  0  or  f  =  f 
y  yy  yy  v 


0  or  exponential 
I  approach 


f  =  0 


fx  =  0 


f  =  -1/2 


f  =  0 
'  xx 


U  f  -  F0(x).  “  fy  =  0 

Figure  5.  Boundary  conditions  for  rapid  expansion  channel 

The  other  boundary  conditions  are  as  in  Fig.  5,  where  the  several  alternatives 
that  we  have  used  for  the  downstream  boundary  condition  are  shown.  If  the 
variable  grid  derivative  coefficients  are  used,  the  table  of  coefficients  for 
the  boundary  conditions  becomes,  for  example; 


S 


•  tr  - 


where  q^  and  r^*  are  given  by  (4.6),  if  we  choose  the  downstream  condition 
f  =  =  0  as  represented  by  (4.4).  This  condition  we  call  the  boundary- 

layer  condition  since  it  implies  that  the  streanwise  second  and  third  deriva¬ 
tives  of  the  streamwise  velocity  are  small  near  the  downstream  boundary.  We 
compare  its  affectiveness  with  that  of  the  following  nonlinear  downstream  con¬ 
dition  which  models  the  exponential  approach  to  the  Poiseuille  profile 

3 

investigated  by  Wilson  . 

If  we  assume  that  the  decay  into  the  fully  developed  Poiseuille  profile 
fp(x)  as  y  -*•  oo  is  exponential,  so  that  we  can  write  for  the  stream  function 
f(x,y) 

f(x,y)  *  fp(x)  +  e'oty4,(x)  (8.1) 

it  can  be  shown  (Wilson  )  that  a  and  4>(x),  which  depend  on  R,  satisfy  an 
equation  similar  to  the  Orr-Sommerfeld  equation,  and  Wilson  has  obtained  an 
approximate  formula  for  a  when  R  is  large.  When  R  is  not  large  enough, 
so  that  we  do  not  know  a  accurately,  we  can  impose  the  asymptotic  condition 

(1)  by  using  nonlinear  boundary  conditions  obtained  by  eliminating  a,  as 
shown  in  the  last  progress  report.  If  y^  is  the  last  y-station  considered 
and  we  write 


the  two  nonlinear  boundary  conditions  consistent  with  (1)  read: 


fN+l  *  fp  +  (fN  VG1  ’  fN  =  fp  +  ^fN-l  VG2 

The  Newtonized  forms  of  these  read 

Vl  ~  a*N  “  b*N-1  =  fp  +  (fN  “  VG1  “  fN-l 

♦n  "  C*N-1  “  d*N-2  =  fp  +  (fN-l  “  fp^G2  ”  fN 

where 

*-  0  ♦  E,)G,,  b  -  "V 

c  ■  0  *  E2)G2,  d  .  -E2G2(fN_,  -  fp)/(fN.2  -  fp) 
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These  nonlinear  boundary  conditions  have  now  been  implemented  in  the 
biharmonic  code  and  compared  with  simpler  downstream  boundary  conditions  for 
the  sudden  expansion  problem.  The  results  in  the  table  below  show  the  decay 
of  the  centerline  velocity  for  a  21  x  31  grid  with  yM„  =  3.0  and  for  the 
same  grid  with  ym_„  =  2.46,  1.94,  1.48  and  1.07,  i.e.  with  the  asymptotic 
conditions  imposed  at  j  =  28,  25,  22  and  19  instead  of  j  =  31 .  For  the 
nonlinear  boundary  conditions  the  fp(x)  distribution  used  was  the  one 
appropriate  to  the  grid,  i.e.  it  was  obtained  by  solving  fxxxx  =  0  on  the 
given  x  grid.  This  is  the  approximation  to  the  fully-developed  profile  that 
one  would  expect  the  computed  profiles  to  tend  to  as  y  «>.  For  the  simpler 
downstream  conditions  we  used  fyy  =  f  =  0. 

The  comparison  shows  that  for  graphical  accuracy  it  is  quite  adequate  to 
impose  the  nonlinear  conditions  at  y  =  1.48  instead  of  y  =  3.0;  thus  comparing 
values  at  x  =  1.6  with  those  obtained  for  xmflx  =3.0,  we  see  that  the  non¬ 
linear  conditions  give  an  error  of  0.0085  whereas  for  the  simpler  condition  the 
error  is  about  -0.046. 

One  may  also  observe  from  these  results  the  rate  at  which  the  effect  of 
the  choice  of  boundary  condition  decays  upstream.  For  example,  if  either 
condition  is  imposed  at  x  =  1.48,  the  difference  could  hardly  be  detected  to 
graphical  accuracy  at  x  =  1.07. 

Runs  so  far  made  on  this  problem  have  been  mainly  with  R  =  41.3  and  have 
started  from  an  initial  guess  consisting  of  a  linear  interpolation  between  the 
inlet  conditions  and  the  fully  developed  downstream  profile.  The  iteration 
strategy  has  normally  been  a  maximum  of  2  Newtons  followed  by  "Chords" 
and  the  monitored  iterative  corrections  printed  out  were,  for  the  finest 
grids  considered  so  far  (21  x  61),  for  example, 

0.19,  0.34  x  10'1,  -0.31  x  10'2,  0.64  x  10"3,  -0.72  x  10"4,  0.79  x  10"5 

Thus  here  2  Newtons  +  4  Chords  were  required  to  achieve  the  required  accuracy 
of  1/2  x  10~4  in  f,  and  this  took  about  0.98  minutes  CPU  time  on  the  IBM  370. 

The  development  of  the  longitudinal  velocity  profiles  from  the  inlet  condi¬ 
tions  to  the  downstream  fully-developed  profile  is  shown  in  Fig.  6  for  the  21  x 
31  grid  In  comparison  with  Durst  et  al.  experimental  data4.  For  the  present  cal¬ 
culation,  the  experimental  data  at  the  inlet  was  served  as  the  initial  profiles, 


which  were  far  from  the  fully  developed  and  of  which  the  corresponding  Reynolds 
number  in  our  notation  was  37.3  instead  of  41.3  for  the  corresponding  fully- 
developed  inlet  conditions.  The  agreement  between  the  calculated  and  experimental 
results  is  very  good;  a  small  discrepancy  of  the  velocity  at  the  centerline  in  the 
inlet  region  could  be  caused  by  the  inaccuracy  in  reading  the  experimental  data, 
particularly  at  the  inlet  which  was  plotted  in  a  small-scale  figure. 


Centerline  Velocities  with  Various  Downstream  Conditions  (R  =  41.3) 


X 

exp 

BL 

exp 

BL 

exp 

BL 

exp 

BL 

0.637 

3.0856 

3.0856 

3.0857 

3.0857 

3.0857 

3.0856 

3.0857 

3.0856 

0.733 

2.9283 

2.9283 

2.9284 

2.9284 

2.9284 

2.9284 

2.9284 

2.9283 

0.837 

2.7667 

2.7667 

2.7667 

2.7667 

2.7668 

2.7668 

2.7668 

2.7665 

0.949 

2.6046 

2.6047 

2.6047 

2.6047 

2.6047 

2.6047 

2.6048 

2.6040 

1.070 

2.4465 

2.4465 

2.4464 

2.4465 

2.4465 

2.4465 

2.4468 

2.4451 

1.198 

2.2963 

2.2964 

2.2963 

2.2963 

2.2964 

2.2963 

2.2973 

2.2913 

1.333 

2.1576 

2.1578 

2.1577 

2.1577 

2.1578 

2.1573 

2.1600 

2.1472 

1.476 

2.0333 

2.0334 

2.0334 

2.0334 

2.0334 

2.0327 

2.0407 

1.9955 

1.626 

1.9248 

1.9249 

1.9249 

1.2950 

1.9252 

1.9211 

1.782 

1.8327 

1.8328 

1.8328 

1.8326 

1.8331 

1.8254 

1.944 

1.7564 

1.7564 

1.7564 

1.7562 

1.7585 

1.7262 

2.112 

1.6946 

1.6947 

1.6947 

1.6927 

2.284 

1.6458 

1.6457 

1.6457 

1.6424 

2.459 

1.6079 

1.6078 

1.6082 

1.5910 

2.638 

1.5790 

1.5781 

2.818 

1.5574 

1.5561 

3.000 

1.5414 

1.5340 
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9.0  TURBULENT  FLOW  APPLICATIONS 


Because  of  the  rapid  variation  near  the  wall,  we  replace  the  no-slip  condi¬ 
tion  at  a  suitably  small  distance  from  the  wall.  The  law  of  the  wall  is  the 
well-known  log  law 


zr=  < 1og  (7Ly)  +  c  o-d 

T 

where  u  =  /T7p  and  k  and  c  are  appropriate  constants.  In  the  law-of- 

T  W 

the  wall  region  50v/ut  <y  <  0. 1 6  we  have 


3U  =  ^  1 

3y  <  y 


or 


u 

T 


(9.2) 


So,  eliminating  u  ,  we  have 

T 

U  =  y  ~  log  eKCy2  1^1  (9.3) 

If  all  variables  are  nondimensionalized  this  becomes 

u  =  y  |p-  log  (a*y2  |^)  (9.4) 

where  a  =  Ke<c  =  3.3  and  y  >  50/R^  where  Rw  =  u^Jt/v.  If  we  define  the 
function  W  by 

W(t)  =  yt  log  (aRy2t)  (9.5) 

the  law  of  the  wall  in  terms  of  the  stream  function  F  reads 

Fy  •  “(Fyy)  (9-6) 

In  terms  of  the  variable  grid  difference  formulation  with  y  =  Y(7)  we  have 

b2f0  +  +  b4f2  =  W(d2fo  +  d3f1  +  d4f2)  (9.7) 

If  we  wish  to  express  fQ  in  terms  of  f.j  and  f2>  as  would  be  required  to 
express  the  boundary  condition  in  the  form  (4.1),  we  use  Newton's  method  again 
and  from  an  approximation  TQ  determine  a  better  one  by 

b,7  +  b,f,  +  b.f,  -  W(d,7n  +  d,f,  +  daf9) 

f  -7  --2.S - U - LI - LL2 - U - LIL  (9.8) 

0  0  b2-d2W'(d270  +  d3f1  +d4f2) 


30 


wKere 


W (t)  =  y[log(ctRy2t)  +  1] 

We  also  need  the  derivatives  3fQ/3f-]  and  3fQ/3f 2 -  These  are  immediately 
found  by  differentiating  (9.7)  to  be 

3fo  b3  “  d3W<  3fo  b4  -  d4W ' 

3f 1  '  b2  -  d2U'  ’  3f2  "  b2  -  d2Wr 

These  provide  the  values  for  a  and  b  in  (4.2). 


(9.9) 


A  sample  of  results  obtained  for  developing  two-dimensional,  plane  channel 

flow  is  presented  on  Fig.  7.  The  Reynolds  number  corresponds  to  experimental 

5 

results  of  Comte-Bellot  and  calculated  results  may  be  compared  with  measure¬ 
ments  at  values  of  x/D  from  20,  where  the  measurements  were  used  as  initial 
conditions  for  the  calculations,  to  59.  The  calculated  results  deviate  from 
measurements,  particularly  in  the  near-wall  region,  as  the  flow  develops  from 
the  initial  condition  and  tends  towards  the  measurements  with  further  downstream 
distance.  The  discrepancies  are  undoubtedly  related  to  the  specification  of 
zero  cross-stream  velocity,  in  conjunction  with  the  law  of  the  wall,  at  the 
initial  station.  This  leads  to  mass  continuity  not  being  initially  satisfied 
and  some  distance  downstream  is  required  before  the  consequences  disappear. 

This  erroneous  assumption  can  be  removed  with  some  effort  but  this  was  not  done 
here  since  the  main  purpose  was  to  develop  and  evaluate  numerical  aspects  of  the 
code. 


31 


.9 

.8 


O 


O 


40.5 


Figure  7 
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.  Computed  turbulent  velocity  profiles  In  a  plane  channel; 
comparison  with  measurements  of  Comte-Bellot  (1961). 


10.0  REFERENCES 


1.  Cebeci ,  T. ,  Hirsh,  R.S.,  Keller,  H.B.  and  Williams,  P.G.:  Studies  of 
Numerical  Methods  for  the  Plane  Navier-Stokes  Equation.  Report  No. 

MDC  J8525,  1979. 

2.  Isaacson,  E.  and  Keller,  H.B.:  Analysis  of  Numerical  Methods.  Wiley,  1966. 

3.  Wilson,  S.:  The  Development  of  Poiseuille  Flow.  JFM,  Vol.  38,  1969, 
pp.  793-806. 

4.  Durst,  F.,  Melling,  A.  and  Whitelaw,  J.H.:  Low  Reynolds  Number  Flow  Over 
a  Plane  Symmetric  Sudden  Expansion.  JFM,  Vol.  64,  1974,  pp.  111-128. 

5.  Comte-Bellot,  G. :  Ecoulement  turbulent  entre  deux  parois  paralleles. 

Publ .  Sci.  et  Tech,  du  Ministere  de  1 'Air,  No.  419,  1965. 


33 


T 

r 


Appendix  A 
B I HARMONIC  ADI 

Conte  (1958)  has  described  an  ADI  technique  for  solving  the  steady 

biharmonic  equation  with  boundary  condition  corresponding  to  the  bending  of 

4 

a  plate  without  clamped  edges.  If  v  $  =  0  is  replaced  by  difference  equa¬ 
tions  in  the  form 

(X  +  Y  +  Z )  4>  =  0 

where  X*  replaces  <}>xxxx>  Y*  replaces  <j>yyyy  and  Z*  replaces  2<(>xxyy 
and  the  boundary  conditions  have  been  used  to  eliminate  exterior  <f>  values, 
then  Conte's  scheme  can  be  written 

4  =  4(k)  -  «k+1(X^  +  Y*(k)  +  Z*(k)) 

*(k+1)  =  *  -ak+1(Y^(k+1)  -  Y*{k)) 

where  «k+^  is  an  iteration  parameter  which  can  be  chosen  to  accelerate 
convergence.  Conte  proves  convergence  and  shows  that  a  good  choice  of  cyclic 
values  for  a  20  x  20  grid  is 

ak  =  (0.2)1-k/16  (k  =  1,2 . 8) 

We  are  interested  in  first  normal  derivative  boundary  conditions 
(clamped  edge)  rather  than  2nd  normal  derivative  conditions  (unclamped  edge), 
but  the  same  scheme  can  be  tried.  However,  no  convergence  results  appear  to 
be  available. 

If  we  drop  the  k  index  and  denote  the  new  41  by  and  also  introduce 
a  nonzero  right-hand  side,  the  above  scheme  can  be  written  in  a  form  suitable 
for  computation  as 

(I  +  <*X)$  ■  —  a((Y  +  Z)$  —  r] 

(I  +  aY)($*  —  4)  =  4  —  $ 

Each  stage  here  requires  sweeps  involving  the  solution  of  a  quin-diagonal 
system.  A  subroutine  QUI SOL  was  written  to  solve  a  general  quin-diagonal 
system  of  the  form 

Vj-2  *  Vj-1  *  cj*j  *  djVi  *  ■  rj  0  * 1 . ") 
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with  a-j  =  -  t>i  =  0,  dn  =  en  =  en_-j  =  0.  This  was  used  to  implement  the 

above  scheme  with  a  view  to  using  it  to  solve  approximately  the  Newton  equa¬ 
tions  of  section  2  with  very  little  storage  requirements.  It  was  tested  on 
a  6x6  system  whose  right-hand  sides  were  generated  as  the  row  sums  of 
the  coefficients,  so  that  the  exact  solution  was  known,  namely  a..  =  1. 

*  J 

Convergence  was  excellent  and  the  exact  solution  obtained,  but  when  the  the 
dimensions  were  increased  to  11  x  11,  convergence  became  very  slow  and 
adjustment  of  the  cyclic  iteration  parameter  did  not  have  much  effect.  The 
kind  of  convergence  that  Conte  predicted  was  not  obtained  at  all,  so  we 

conclude  that  the  chanqe  of  boundary  condition  makes  a  vital  difference, 
and  that  some  modification  to  the  scheme  would  certainly  be  needed. 


i 
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APPENDIX  B 
BOX  ADI  SCHEME 

The  Box  ADI  scheme  was  originally  proposed  for  dealing  with  the  coupled 
elliptic  equations  which  occur  In  comer  boundary  layers.  Since  these  reduce 
essentially  to  the  two-dimensional  boundary- layer  equations  at  large  distances 
from  the  comer  the  view  was  taken  that  it  would  be  generally  beneficial  to 
construct  a  scheme  that  would  reduce  to  the  familiar  two-dimensional  Box  scheme 
far  away  from  the  comer.  A  fully  three-dimensional  Box  scheme  has  several 
disadvantages,  not  least  of  which  is  the  proliferation  of  dependent  variables, 
so  we  proceed  as  though  a  standard  ADI  scheme  based  on  central  differences  was 
intended  then  use  the  ordinary  Box  scheme  for  the  boundary  value  problems  that 
have  to  be  solved  on  each  grid  line.  This  means  that  slightly  different  differ¬ 
ence  approximations  are  being  used  in  each  direction,  but  they  are  all  second 
order  so  this  should  not  matter  too  much;  in  fact  there  is  one  way  in  which  it 
may  be  an  advantage,  which  is  best  explained  after  the  details  of  the  scheme 
have  been  described. 

We  considered  first  a  simple  example  to  test  the  scheme,  namely 

Uxx  +  uyy  +  aux  +  buy  +  cu  =  f(x»y)  +  aut 

where  f(x,y)  was  chosen  so  that  the  equation  had  a  given  known  solution  and 
appropriate  boundary  values  were  imposed  as  a  rectangle.  A  grid  which  may  be 
nonuniform  in  each  direction,  is  set  up  to  cover  the  rectangle  as  in  Fig.  3. 

In  the  standard  ADI  scheme  we  first  consider  sweeps  In  the  x  direction  and 
determine  u  at  t  +  %at  to  satisfy 

“xx  ♦  »“x  +  <?  C  -  fjr>“  “  f(x"*>  -  “yy  -  byy  -  c  *  ZT>“  <B1  > 

on  each  x  grid  line  j  *  2,  ...»  n-1,  where  u  is  the  known  field  at^the 
previous  time  step  t.  We  now  sweep  in  the  y-direction  to  solve  for  u 

“yy  ft>“  *  f<x^>  -  “xx (K) 

on  each  y-grid  line  1  *  2,  ....  n-1.  We  then  take  u  as  the  updated  field 
for  t  +  at.  Since  we  are  Interested  In  the  stea<1y  state,  we  repeat  these 
sweeps  until  convergence. 
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In  the  standard  ADI  scheme  we  solve  the  boundary  value  problems  associated 
with  (Bl)  and  (B2)  by  central  differences  and  a  tridiagonal  solver,  with 
the  right-hand  side  derivatives  evaluated  by  central  differences.  In  the  Box 
ADI  scheme  we  solve  them  by  the  Box  method  and  a  S0LV2  subroutine,  but  still 
evaluate  the  right-hand  side  derivatives  by  central  differences.  For  (Bl),  for 


example,  we  introduce  v  =  ux  and  consider  the  first-order  system 


ux  -  v  =  0 


(B3) 


vx  +  av  +  (|  -  ||)u  =  f(x,y)  -  u**  -  bu*  -  (f  +  ||)u 


where  u**,  u 


2  At 
*  are  u 


ux  evaluated  by  central  differences. 


If  h..  =  xi+i  —  x^,  then  the  Box  scheme  approximations  yield  for  i  =  2, 
•  %  m 


(ui  ”ui-i>  "  7  (vi  +  vi-l)  = 


l 


-h^7  <Vi  -  V1-1J  +  f  (VT  +  Vl}  +  -  zf)(G1  +  Ui-1}  =  ri 


(B4) 


where 


r,  -  $  [f(x,_y)  *  -fa'  ♦  -t*! )  -  |(uf  *  -\ (f  *  |f)  CUl  ♦  V)) 


These  equations,  together  with  the  boundary  conditions,  are  solved  either  by  a 

a 

block  tri-diagonal  solver  or  a  band  solver.  We  can  overwrite  the  u  on  the  u 
provided  we  have  stored  temporarily  the  old  u's  needed  to  calculate  the  u**, 
u*  for  the  next  grid  line. 


Having  performed  those  calculations  for  each  x  grid  line,  i.e.,  j  =  2, 
....  n-1,^  we  change  to  sweeping  in  the  y-dlrectlon.  We  can  use  v  again,  but 
now  for  Uy.  The  system  we  now  solve  along  each  x  grid  line  is  the  same  as 
(B3)  except  that  the  subscript  x  becomes  subscript  y,  a  and  b  are  inter¬ 
changed  and  u**,  U*  now  mean  derivatives  in  the  x*direction  evaluated  by 
central  differences.  Further  we  interpret  0  as  0  and  u  and  u,  which 

is  natural  anyway  since  u  will  have  been  overwritten.  The  Box  equations  for 
the  y-dlrectlon  follow  by  similar  interpretations  and  also  changing  1  to  j 


and  Interpreting  hj  as  y^-y^. 
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One  might  expect  that,  once  the  iterations  had  converged  to  a  level  cor¬ 
responding  to  the  differences  between  the  two  approximation  schemes,  the  iter¬ 
ative  changes  would  not  decrease  further  and  that  this  level  would  therefore 
give  an  indication  of  the  truncation  error;  thus  we  would  have  an  automatic 
criterion  for  terminating  the  ADI  iterations.  However,  on  the  simple  examples 
tested  the  iteration  seemed  to  converge  to  normal  working  accuracy,  approximately 
six  significant  figures. 

When  the  scheme  was  tried  on  the  vorticity/stream-function  formulation  of 
the  Navier-Stokes  equations,  the  situation  was  not  so  satisfactory.  We  now 
have  two  coupled  equations,  of  course,  and  moreover  we  have  two  boundary  condi¬ 
tions  on  one  variable  but  none  on  the  other,  so  the  situation  is  somewhat 
different.  However,  it  was  hoped  that  by  solving  the  two  equations  simultaneously 
on  each  grid  line  by  the  Box  method,  the  boundary  condition  problem  would  be 
automatically  dealt  with. 


If  f  is  the  stream  function  and  q  the  vorticity,  the  equations  read 

f  +  f  -  q  +  af . 
xx  yy  H  t 

i  (q  +q  )  =  f  q  — fq  +  q. 

R  VHxx  Hyy'  yHx  xHy  Ht 


(B5) 


where  aft  is  a  fictitious  time-dependent  term  introduced  to  help  convergence, 
if  necessary. 

A  A 

For  the  x-sweep  equations  we  Introduce  g  =  fx  and  s  =  qx»  and  the  equa¬ 
tions  corresponding  to  (B3)  then  read 

f„  -  g  •  0 


(B6) 


qx“s 


sx-(RVS“f  q+  (Rq*)g  =  ~  q  -  q** 


These  are  discretized  by  the  Box  scheme  and  solved  by  a  block  tri-diagonal 
solver  or  a  band  solver  In  exactly  the  same  way  as  before,  and  similarly  for 
the  y-sweep  equations.  Note  that  we  do  not  require  Initial  guesses  for  g 
and  s  because  they  do  not  occur  as  coefficients  or  in  the  right-hand  side. 
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Successive  pairs  of  alternate  sweeps  using  this  scheme  should  gradually 
update  the  f  and  q  fields  to  the  solutions  of  the  steady-state  equations. 

In  fact,  this  occurred,  when  the  scheme  was  tried  on  the  simple  Poiseuille 
flow  with  the  initial  guess  taken  to  be  the  exact  solution  plus  a  perturbation 
which  was  zero  at  the  comers. 

Now  there  is  a  slight  difficulty  at  the  comers.  In  the  x-sweeps  we 
update  the  vorticity  on  the  whole  of  each  grid  line  for  j  =  2,  ....  n-1. 

Thus  we  obtain  updated  q  values  along  the  left-  and  right-hand  boundaries, 
but  not  at  the  comers.  When  we  proceed  to  the  y-sweeps,  we  update  q  along 
the  top  and  bottom  boundaries,  but  not  at  the  comers.  Thus  some  other  means 
must  be  found  for  updating  the  corner  vorticity  values  since  they  are  needed  to 
evaluate,  for  example,  the  q**  boundary  values.  We  chose  to  update  the  value 
of  q  at  a  comer  by  taking  the  average  of  the  two  values  obtained  by  extrap¬ 
olating  quadratically  along  each  of  the  adjacent  sides.  With  this  procedure 
the  Iterations  converged  for  the  simple  Poiseuille  flow,  even  when  the  pertur¬ 
bation  from  the  exact  solution  was  not  zero  at  the  comers.  It  was  found  that 
the  fictitious  time  derivative  term  in  (B5)  was  not  needed,  so  a  was  set  to 
zero. 


When  the  scheme  was  applied  to  the  channel  inlet  problem  with  a  uniform 
parallel  entry  velocity,  no  convergence  could  be  obtained  at  all.  Here  the 
initial  guess  chosen  was  fully  developed  flow  everywhere  except  right  at  the 
Inlet,  where  the  entry  conditions  were  Imposed.  Because  there  are  infinite 
vorticity  values  at  the  comers  In  this  case.  It  seems  highly  likely  that  It  is 
the  comer  problem  that  Is  affecting  the  rest  of  the  solution.  We  hope  to  try 
the  scheme  on  a  more  realistic  problem  with  Infinite  vorticity  values  as  a  part 
of  the  future  program. 
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